home *** CD-ROM | disk | FTP | other *** search
/ Resource for Source: C/C++ / Resource for Source - C-C++.iso / codelib7 / v_09_09 / 9n09012a < prev    next >
Encoding:
Text File  |  1995-11-01  |  1.9 KB  |  73 lines
  1. /* _Sin function */
  2. #include "xmath.h"
  3.  
  4. /* coefficients */
  5. static const double c[8] = {
  6.     -0.000000000011470879,
  7.      0.000000002087712071,
  8.     -0.000000275573192202,
  9.      0.000024801587292937,
  10.     -0.001388888888888893,
  11.      0.041666666666667325,
  12.     -0.500000000000000000,
  13.      1.0};
  14. static const double s[8] = {
  15.     -0.000000000000764723,
  16.      0.000000000160592578,
  17.     -0.000000025052108383,
  18.      0.000002755731921890,
  19.     -0.000198412698412699,
  20.      0.008333333333333372,
  21.     -0.166666666666666667,
  22.      1.0};
  23. static const double c1 = {3294198.0 / 2097152.0};
  24. static const double c2 = {3.139164786504813217e-7};
  25. static const double twobypi =
  26.     {0.63661977236758134308};
  27. static const double twopi =
  28.     {6.28318530717958647693};
  29.  
  30. double _Sin(double x, unsigned int qoff)
  31.     {   /* compute sin(x) or cos(x) */
  32.     switch (_Dtest(&x))
  33.         {
  34.     case NAN:
  35.         errno = EDOM;
  36.         return (x);
  37.     case 0:
  38.         return (qoff ? 1.0 : 0.0);
  39.     case INF:
  40.         errno = EDOM;
  41.         return (_Nan._D);
  42.     default:    /* finite */
  43.          {  /* compute sin/cos */
  44.         double g;
  45.         long quad;
  46.  
  47.         if (x < -HUGE_RAD || HUGE_RAD < x)
  48.             {   /* x huge, sauve qui peut */
  49.             g = x / twopi;
  50.             _Dint(&g, 0);
  51.             x -= g * twopi;
  52.             }
  53.         g = x * twobypi;
  54.         quad = (long)(0 < g ? g + 0.5 : g - 0.5);
  55.         qoff += (unsigned long)quad & 0x3;
  56.         g = (double)quad;
  57.         g = (x - g * c1) - g * c2;
  58.         if ((g < 0.0 ? -g : g) < _Rteps._D)
  59.             {   /* sin(tiny)==tiny, cos(tiny)==1 */
  60.             if (qoff & 0x1)
  61.                 g = 1.0;    /* cos(tiny) */
  62.             }
  63.         else if (qoff & 0x1)
  64.             g = _Poly(g * g, c, 7);
  65.         else
  66.             g *= _Poly(g * g, s, 7);
  67.         return (qoff & 0x2 ? -g : g);
  68.          }
  69.         }
  70.     }
  71. /* End of File
  72.  
  73.